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1.  Introduction 

Amorphous  silicon  dioxide  (a-Si02)  comprises  the  gate  and  passivation  oxide  layers  in  most  MOS 
devices.  However,  the  current  understanding  of  defects  induced  in  a-Si02  due  to  high-energy  radiation 
in  space  is  inadequate  for  making  predictions  about  various  options  for  device  fabrication.  It  is  also 
inadequate  for  interpretation  of  results  of  radiation  testing  of  critical  components  of  Air  Force  Space 
Systems.  Theoretical  modelling  plays  an  important  role  in  developing  models  and  elucidating 
properties  of  these  defects.  However,  computational  methods  currently  available  for  ab  initio 
calculations  of  defects  in  amorphous  solids  are  not  fully  adequate  for  predictive  modelling.  They  fall 
into  two  major  groups.  In  the  first  group,  an  infinite  amorphous  structure  is  represented  by  periodic 
translation  of  an  amorphous  unit  cell.  These  methods  often  employ  ab  initio  Molecular  Dynamics  to 
generate  an  amorphous  structure  and  then  use  the  same  computational  technique  in  static  mode  to 
study  defects  in  this  structure  (see,  for  example,  refs.  [1-3]).  In  the  second  group,  amorphous  structure 
is  represented  by  a  finite  cluster  saturated  at  the  border  by  Hydrogen  atoms  or  model  pseudo-atoms 
[4],  Small  molecular  clusters  used  in  most  ab  initio  calculations  cannot  account  for  both  the  Madelung 
potential  in  amorphous  structure  and  the  long-range  structural  relaxation  induced  by  a  defect.  These 
clusters  are  either  cut  out  from  a  larger  fragment  of  an  amorphous  structure  (see  section  2)  or 
artificially  constructed  from  general  considerations  to  model  the  defect  enviromnent.  The  periodic 
approach  is  more  consistent  in  a  sense  that  the  same  ab  initio  method  can  be  used  to  produce  the 
amorphous  structure  and  to  study  defect  properties.  However,  a b  initio  calculations  of  periodic 
systems  are  time-consuming  and  therefore  consider  small  periodic  cells.  This  does  not  allow  one  to 
fully  account  for  the  intermediate  and  long-range  order  in  amorphous  structures  [5,6]. 

An  alternative  approach  is  an  embedded  cluster  model  where  single  defect  is  treated  in  infinite 
polarizable  lattice.  A  hybrid  method  developed  in  our  group  at  UCL  [7-11]  combines  quantum 
mechanical  treatment  of  atoms  surrounding  a  defect  with  the  classical  representation  of  the  rest  of  the 
solid  in  a  shell  model  [12].  It  has  been  widely  used  to  study  point  defects  in  ionic  crystals.  This 
approach  naturally  account  for  both  the  Madelung  potential  at  the  site  of  interest  and  the  ionic  and 
electronic  polarisation  of  the  system  by  the  defect.  In  the  EOARD  funded  contract  No  F6 1 775-00- 
WE001  we  extended  the  embedded  cluster  method  to  include  ionic-covalent  systems,  such  as 
crystalline  Si02.  Tire  main  aim  of  this  project  was  to  develop  this  method  further  that  it  could  treat 
disordered  (amorphous)  Si02.  This  involved:  1)  developing  of  a  computer  code  for  comprehensive 
analysis  of  amorphous  structures;  and  2)  extending  the  embedded  cluster  approach  to  calculation  of 
localised  defects  in  amorphous  Si02,  and  making  test  applications  to  charged  Hydrogen  defects: 
proton  (H+)  and  H  centers. 

2.  Details  of  calculations 

Unlike  in  crystalline  structures,  all  sites  in  amorphous  structures  are  different  and  one  should  generally 
study  distributions  of  defect  properties  over  many  different  sites.  On  the  other  hand,  some  sites  may 


DiS™BUJ'°N  STATEMENT  A 

ApPr°)edf°r  Public  Release 
Distribution  Unlimited 


4Q  roa-  IZ'3351 


have  specific  properties,  which  make  them  particularly  important  for  applications.  For  example,  three 
and  four-member  rings  in  amorphous  Si02  are  thought  to  be  responsible  for  the  low-energy  optical 
absorption  tail  which  can  be  selectively  excited  by  7.9  eV  photons  of  F2  laser  [13].  Therefore  the 
general  strategy  of  embedded  cluster  calculations  of  defects  in  a-Si02  includes:  i)  generating 
amorphous  structures;  ii)  characterizing  obtained  structures  in  terms  of  distributions  of  various 
structural  parameters,  and  finding  particular  sites  interesting  for  defect  studies;  iii)  constructing 
quantum  clusters  at  different  sites  and  calculating  defect  properties;  iv)  building  a  statistics  of  defect 
properties  over  many  sites  and  producing  distributions  of  defect  properties  (see,  for  example,  ref. 
[14]).  The  main  objective  of  this  project  was  to  build  computational  methods  for  each  stage  of  this 
program  and  to  test  them  on  a  particular  system.  These  methods  are  discussed  in  detail  below. 

Preparation  of  amorphous  structures.  Several  methods  have  been  used  in  the  literature  for 
generating  a-Si02  structures  (see,  for  example,  refs.  [5,6,15]).  In  this  work  we  resorted  to  the  most 
commonly  employed  classical  molecular  dynamics  technique  [16,17].  This  technique  can  be  used  to 
model  melting  and  cooling  process  of  some  crystalline  Si02  polymorph  and  in  this  way  to  mimic  a 
process  of  glass  formation  in  laboratory  conditions.  It  is  particularly  convenient  for  our  applications  as 
it  is  compatible  with  the  embedded  cluster  method  where  the  same  classical  atomistic  simulation 
technique  is  used  for  representing  the  rest  of  the  amorphous  structure  surrounding  a  defect  treated 
quantum-me  chani  c  al  ly . 

The  molecular  dynamics  (MD)  simulations  of  the  glass  formation  were  made  in  the  periodic  model 
using  the  DL_POLY  computer  code  [18].  The  crystalline  P-crystobalite  is  a  convenient  choice  of 
initial  structure  due  to  its  cubic  lattice  cell.  The  classical  inter-atomic  potentials  used  in  this  study  were 
developed  in  ref.  [19]  using  the  results  of  ab  initio  calculations  of  small  Si02  clusters.  These  potentials 
reproduce  accurately  the  relative  stability  of  several  Si02  polymorphs  and  have  been  used  in  several 
previous  studies  of  Si02  glasses  [17].  The  glass  preparation  procedure  is  summarized  below. 

i)  Two  super-cells  of  P-crystobalite  of  192  and  648  atoms  were  heated  relatively  fast  (at 
approximately  4T013  K/s)  at  constant  normal  pressure  to  high  temperature  (Thigh  =  7000K). 

ii)  The  system  was  kept  at  Thigh  for  approximately  100  ps,  which  is  significantly  longer  than  is 
required  to  achieve  thermal  equilibrium. 

iii)  The  equilibrated  system  was  cooled  down  to  OK  or  300K  with  a  relatively  small  cooling  rate, 
which  is  sufficient  to  provide  convergence  of  characteristic  glass  properties  [17].  The  slowest 
cooling  rate  used  in  this  study  was  1012  K/s. 

iv)  Finally,  the  thus  produced  structures  were  relaxed  using  a  static  energy  minimization 
technique. 

Several  amorphous  structures  were  obtained  using  this  procedure  and  different  cooling  rates  and  other 
parameters  of  molecular  dynamics.  The  main  characteristics  of  these  structures  were  analyzed  using 
the  method  described  below  and  are  discussed  in  section  3. 

Analysis  of  amorphous  structures  is  required  to:  i)  compare  different  amorphous  structures;  ii) 
identify  potential  traps  for  electrons,  holes  and  excitons  in  the  ‘non-defective’  amorphous  structure; 
and  iii)  correlate  local  properties  of  amorphous  structure  with  defect  properties,  such  as  optical 
excitation  and  defect  formation  energies.  Out  of  several  amorphous  structures  generated  by  the  method 
described  above  we  have  chosen  only  those  structures  which  had  ideal  connectivity  and  thus 
corresponded  to  the  continuous  random  network  model  of  ideal  glass  [5,6].  The  latter  assumes  the 
ideal  4-2  coordination,  i.e.  each  Si  atom  is  coordinated  by  four  Oxygen  atoms  and  each  O  atom  is 
coordinated  by  two  Silicon  atoms.  The  experimental  interference  and  correlation  functions  deduced 
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from  neutron  diffraction  data  provide  information  about  three  different  ranges  of  order  in  the 
amorphous  structure  [6],  Range  I  (short-range  of  ca.  1-2  A)  describes  the  local  geometrical  parameters 
such  as  Si-0  inter-atomic  distances  and  O-Si-O  tetrahedral  angles.  Range  II  (medium  range  of  ca.  2- 
5  A)  includes  parameters  such  as  Si-Si  inter-atomic  distances  and  torsion  angles,  which  describe 
relative  orientation  of  neighboring  Si04  tetrahedral  units.  Range  III  corresponds  to  long-range  order  at 
ca.  5-15  A  from  the  center  of  interest.  We  should  note  that  both  the  experimental  and  theoretical 
methods  rely  on  a  number  of  approximations.  In  particular,  inter-atomic  potentials  are  fitted  to 
reproduce  crystal  structure  and  relative  energies  of  different  Si02  polymorphs.  These  potentials  are 
best  defined  in  a  region  close  to  the  equilibrium  separation  between  atoms  in  a  crystal  and,  therefore, 
may  not  necessarily  reproduce  the  Si02  melt.  Also  the  typical  size  of  a  system  and  timescale  of 
presently  available  MD  simulations  are  many  orders  of  magnitude  smaller  than  those  in  real  life 
processes,  which  represents  a  principal  difficulty.  Nevertheless,  it  is  presently  accepted  that  thus 
obtained  amorphous  structures  are  representative  of  real  amorphous  systems  and  therefore  it  is 
essential  to  have  a  tool  for  comprehensive  analysis  of  their  properties. 

In  this  project,  we  developed  a  computer  code  PAS  (Portrait  of  an  Amorphous  System)  for  geometric 
and  correlation  analysis  of  crystalline  and  amorphous  Si02  structures.  The  code  is  currently  applicable 
to  systems  with  the  ideal  4-2  coordination  and  can  be  easily  extended  to  systems  with  non-ideal 
coordination  too.  It  inputs  parameters  of  the  unit  cell,  fractional  or  Cartesian  coordinates  of  unit  cell 
atoms  and  values  of  on-site  electrostatic  potential.  It  first  determines  a  table  of  nearest  neighbors 
(TNN)  using  the  “smallest  distance”  criterion  (e.g.  O  atoms  at  the  distance  less  than  1.8  A  from  a  Si 
atom  are  taken  as  its  nearest  neighbors).  The  number  of  characteristics  is  then  calculated  using  the 
TNN: 

i)  distributions  of  Si-0  and  Si-Si  inter-atomic  distances; 

ii)  distributions  of  Si-O-Si  and  O-Si-O  angles; 

iii)  electrostatic  potential  distribution; 

iv)  pair-wise  distribution  functions  g(Si-O),  g(O-O),  g(Si-Si); 

v)  number  of  different  rings  in  the  system. 

The  PAS  code  also  perfonns  analysis  of  the  ring  structure  and  correlates  the  properties  i)  -  iii)  of  ring 
atoms  with  the  size  of  the  ring.  In  addition  the  PAS  code  generates  an  input  for  embedded  cluster 
calculations  using  the  GUESS  code. 

Analysis  of  rings  in  a-Si02  is  often  used  as  a  tool  for  characterizing  the  properties  of  amorphous 
structures.  Structure  of  many  ion-covalent  and  covalent  materials,  in  general,  and  Si02  polymorphs,  in 
particular,  allows  for  closed  paths  (rings)  within  the  system.  The  number  and  size  of  rings  in  different 
Si02  polymorphs  is  associated  with  their  structure  and  density  [5].  Statistics  and  structure  of  rings  in  a- 
Si02  is  of  considerable  interest  because  small  planar  rings  are  thought  to  determine  the  mobility  edge 
and  the  optical  absorption  tail  in  real  glasses  [13,20]  and  larger-size  rings  are  responsible  for  gaseous 
permeability  and  diffusion  of  defect  species.  The  small  3-  and  4-member  rings  have  also  been 
associated  with  specific  peaks  (Dt  and  D2)  in  the  Raman  spectrum  [2,3]. 

There  are  two  most  commonly  used  definitions  of  a  ring.  One  of  them  defines  a  ring  associated  with 
an  atom  Si*  as  a  shortest  path  Si;  -  Oy  - Oik  -  Sii,  where  Ojj  and  0*  are  nearest  neighbors  of  Sij  and  j 
*  k.  According  to  this  definition  each  Si  atom  has  6  rings  associated  with  it.  We  have  adopted  another 
definition  in  which  a  ring  is  defined  as  a  closed  path  without  shortcuts,  i.e.  without  any  path  from  one 
atom  of  the  ring  to  the  other  being  shorter  than  the  ring  path.  The  two  definitions  are  not  identical,  but 
the  second  definition  seems  more  meaningful  because  it  regards  a  ring  as  an  entity,  a  particular  feature 
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of  the  glass  structure,  while  the  first  one  considers  rings  in  the  mathematical  sense  as  a  number  of 
closed  paths  associated  with  each  atom. 


Fig.  1.  Tree  of  neighbors  for  a  root  Si  atom  at  the  top  of  the  figure.  Black  and  white  circles  arranged  in  levels  are 
Si  and  O  atoms  respectively.  The  level  number  is  shown  on  the  left;  values  on  the  right  show  the  total  number  of 
atoms  in  each  level.  The  bold  line  shows  a  closed  3-member  path  found  for  the  root  Si  atom:  the  path  includes  Si 
atom  at  level  6  (the  same  atom  as  at  level  0)  and  two  other  Si  atoms  at  levels  4  and  2  and  corresponding  Oxygen 
atoms. 

Analysis  of  the  ring  statistics  in  a-Si02  is  a  complex  computational  problem.  To  explain  the  algorithm 
of  ring  search  implemented  in  the  PAS  code  we  will  refer  to  Figs.  1  and  2.  The  search  for  rings  begins 
with  construction  of  a  tree  of  neighbors  for  each  Si  atom,  as  shown  in  Fig.  1.  In  such  a  tree  the  ‘root” 
Si  atom  occupies  level  0  and  further  levels  are  occupied  by  oxygen  and  silicon  atoms  in  sequential 
order.  The  tree  is  built  so  that  to  include  sufficiently  large  number  of  shells  of  neighbors  (up  to  20 
levels  in  our  calculations).  Each  atom  in  each  Si  level  (SiieVei,vemx)  is  then  compared  with  the  Si  atom  at 
the  root  of  the  tree.  A  closed  path  via  the  root  atom  is  found  if  atom  Siievei,vertix  atom  has  the  same 
fractional  coordinates  as  the  root  atom.  An  example  of  such  a  path  is  shown  in  Fig.  1.  The  procedure  is 
repeated  for  each  Si  atom  in  the  unit  cell.  All  found  cycles  are  accumulated  for  further  analysis. 


Fig.  2.  Closed  paths  removed  from  the  list  of  candidates  for  rings.  (Only  Si  atoms  are  shown.)  a)  Paths  a  and  b 
differ  by  the  direction  of  cycling  only,  one  of  them  is  removed,  b)  At  least  one  atom  appears  twice  in  the  same 
path,  c)  The  path  connects  the  same  atom  in  different  (not  necessarily  adjacent)  super-cells,  d)  Paths  c  is  shorter 
than  paths  a  or  b  (or  both). 

At  the  next  step  the  structure  and  topology  of  all  found  cycles  are  analyzed.  By  construction,  the  list  of 
all  cycles  found  at  the  previous  step  is  excessive,  i.e.  it  contains  many  equivalent  cycles  and  those, 
which  do  not  fall  under  the  definition  of  a  ring.  To  exclude  equivalent  paths  we  note  that  both 
clockwise  and  anti-clockwise  paths  will  be  found  for  each  cycle  (paths  a  and  b  in  Fig.  2(a)).  In 
addition,  each  cycle  of  length  n  (w-member  cycle,  where  n  is  the  number  of  Si  atoms)  will  be  found  n 
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times.  This  is  because  each  w-member  cycle  will  be  found  in  n  different  trees  associated  with  each  Si 
atom  in  the  cycle.  Among  remaining  cycles  the  following  ones  have  to  be  excluded:  i)  those  which 
have  any  atom  appearing  more  than  once  (as  in  Fig.  2(b));  ii)  those  which  connect  the  same  Si  atom  in 
different  unit  cells  (as  in  Fig.  2(c));  and  iii)  those  which  have  shortcuts  (path  c  in  Fig.  2(d)  is  shorter 
than  paths  a  and  b).  The  remaining  cycles  are  called  rings. 

The  embedded  cluster  method  developed  in  our  group  and  implemented  in  the  GUESS  computer 
code  allows  us  to  study  point  defects  in  crystalline  and  amorphous  materials  combining  quantum 
mechanical  treatment  of  a  defect  and  its  nearest  surrounding  with  the  classical  shell  model  treatment  of 
the  rest  of  the  solid.  The  GUESS  method  and  details  of  its  applications  to  point  defects  in  crystalline 
Si02  were  described  in  the  final  report  to  the  contract  NoF61775-00-WE001  and  in  refs.  [9-11].  Here 
we  highlight  only  the  features  specific  to  the  present  studies. 

In  the  GUESS  method,  the  infinite  system  with  a  single  point  defect  is  divided  into  several  regions.  A 
spherical  region  I  includes:  i)  a  quantum  cluster  with  a  defect  and  surrounding  atoms  treated  quantum 
mechanically;  ii)  an  interface  region,  which  connects  the  quantum  cluster  with  the  rest  of  the  solid 
treated  classically;  iii)  a  classical  region,  which  includes  up  to  several  hundred  atoms.  Region  I  is 
surrounded  by  a  finite  region  II,  which  is  treated  atomistically,  and  region  III,  which  is  treated  in  the 
approximation  of  polarisable  continuum.  Region  III  conforms  geometrically  to  the  boundary  between 
regions  I  and  II  but  extends  to  infinity.  The  classical  ions  in  regions  I  and  II  are  treated  in  the  shell 
model  and  interact  between  themselves  via  inter-atomic  potentials  [19].  Note  that  the  same  inter¬ 
atomic  potentials  are  used  to  generate  amorphous  structures  using  classical  MD  simulations.  Both 
quantum  and  classical  atoms  in  region  I  are  allowed  to  relax  in  the  course  of  calculations.  Atoms  in 
region  II  are  kept  fixed  in  their  ideal,  bulk  positions  and  provide  correct  variations  of  the  electrostatic 
potential  inside  the  region  I.  Region  III  is  used  to  calculate  the  polarization  energy  of  the  infinite 
lattice  due  to  the  presence  of  a  defect  in  region  I  using  the  Mott-Littleton  approach.  One  of  the  main 
advantages  of  this  approach  is  that  it  allows  us  to  calculate  forces  on  quantum  mechanical  and 
classical  ions  and  simultaneously  optimize  their  positions  using  an  effective  energy  minimization 
scheme. 

Embedded  cluster  calculations  of  amorphous  silica  differ  significantly  from  these  for  crystalline 
material.  First,  we  note  that  all  atoms  in  amorphous  structures,  in  principle,  have  different  local 
structure  of  their  environment.  Therefore,  it  is  often  impossible  to  use  the  same  QM  cluster  to  describe 
properties  of  the  same  defect  at  different  sites.  Then,  comparison  of  results  for  the  same  or  different 
defects  can  only  be  feasible  if  the  dependence  of  the  defect  properties  on  the  cluster  size  and  shape  is 
less  significant  than  their  dependence  on  the  atomic  site.  In  order  to  estimate  an  error  introduced  by 
the  finite  size  of  the  QM  cluster,  and  to  establish  grounds  for  comparison  of  the  results  obtained  for 
different  sites  in  an  amorphous  structure,  we  have  studied  the  dependence  of  the  results  on  the  cluster 
size.  To  achieve  that  we  performed  calculations  of  the  H+  impurity  center  in  a-quartz  using  two 
different  clusters:  the  smaller  cluster  Si207Si*6  and  the  larger  cluster  Si5016Si*12  shown  in  Figs.  3(a) 
and  3(b),  respectively.  The  Si207Si*6  cluster  is  too  small  to  calculate  the  H“  center  and  therefore  only 
the  Si5016Si*12  cluster  was  used  for  that. 

Region  I+II  in  embedded  cluster  calculations  may  have  different  shapes  but  should  be  neutral  with 
zero  dipole  moment.  For  calculations  of  H  centers  in  a-Si02  we  have  chosen  a  spherical  region  I+II  of 
the  radius  30.0  A.  This  region  was  constructed  from  the  648  atom  supercell  obtained  in  MD 
simulations  and  was  split  into  Si(01/2)4  tetrahedral  units.  The  latter,  due  to  their  rigid  character,  provide 
a  convenient  choice  of  structural  elements  in  both  crystalline  and  amorphous  Si02.  Two  different 
regions  I+II  were  built  for  a-Si02:  one  was  centered  on  a  3-member  ring  and  the  other  -  on  a  6- 
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member  ring.  Owing  to  the  disordered  structure  the  two  regions  have  different  number  of  atoms  in 
them:  7912  for  the  3-  and  7876  for  the  6-member  ring.  For  the  H  centers  in  alpha-quartz  we  took 
region  I+II  built  from  asymmetric  9-atom  unit  cell  repeated  14  times  along  each  of  its  translation 
vectors.  This  region  contains  almost  25000  atoms.  In  all  cases  a  spherical  region  I  of  the  radius  of  13.0 
A  was  used;  however,  the  number  of  atoms  in  it  was  different  for  different  systems:  659  and  658  for 
the  3-  and  6-member  rings  respectively  and  718  for  alpha-quartz. 

To  describe  accurately  the  structure  and  properties  of  a  defect  associated  with  a  ring,  it  is  convenient 
to  build  the  corresponding  QM  cluster  so  that  it  includes  all  ring  atoms  and  their  nearest  neighbors. 
For  an  ^-member  ring  such  cluster  includes  n  Silicon  atoms  and  n  Oxygen  atoms  which  belong  to  the 
ring,  2 n  O  atoms  which  formally  do  not  belong  to  the  ring  but  coordinate  the  rings’  Si  atoms,  and  2 n 
pseudo-Si  atoms  which  link  the  QM  cluster  with  its  classical  environment.  Such  a  Si,iO„+2„Si*2n  cluster 
is  the  minimal  QM  cluster  required  to  describe  the  ring  in  full.  A  Si309Si*6  cluster  for  a  3-member  ring 
considered  in  this  study  is  shown  in  Fig.  3(c).  A  similar  but  larger  QM  cluster  Si60|8Si*i2  was  used  to 
treat  a  6-member  ring. 


Fig.  3.  (a,  b)  Quantum  mechanical  clusters  Si207Si*6  and  Si50i6Si*i2  used  in  the  calculations  of  H  centers  in  ex- 
quartz  and  a-Si02;  c)  Si309Si*6  -  quantum  cluster  used  to  calculate  the  3-member  ring. 

Not  only  QM  clusters  but  also  location  of  H  centers  within  a  cluster  can  be  different.  Indeed,  in  the 
case  of  a-quartz  all  sites  are  the  same  both  in  terms  of  their  topology  and  environment.  Therefore  it  is 
sufficient  to  consider  H+  and  FT  ions  bonded  to  one  Oxygen  or  Silicon  lattice  site,  respectively.  In  the 
case  of  a  ring  cluster,  however,  all  Si  or  O  sites  of  the  ring  are  topologically  equivalent  but  their 
environment  is  different.  The  latter  should  determine  the  difference  is  properties  of  the  H  centers.  In 
the  present  study  we  did  not  investigate  such  dependence  in  details  but  considered  only  one  site  for 
each  H  center  in  one  3-member  and  one  6-member  ring. 

A  standard  all-electron  6-3 1G  basis  set  was  used  to  describe  lattice  Si  and  O  ions  in  these  calculations. 
The  cc-pVDZ  basis  set  from  the  Gaussian98  library  was  used  for  H  ions.  The  interface  Si*  atoms, 
described  in  details  in  the  final  report  to  the  contract  NoF61775-00-WE001,  were  treated  using  a  local 
effective  core  pseudo-potential  and  a  contracted  sp  basis  set. 

The  incorporation  energies  for  the  Hydrogen  defects  are  calculated  using: 

Emc(Hq)  -  (  E(Hq)  +  E(lattice) )  -  E(lattice  +  Hq), 

where  Hq  is  an  H+  or  H“  ion.  The  electrostatic  potential  inside  the  region  I+II  is  shifted  with  respect  to 
that  obtained  by  Ewald  summation  by  Vshm.  Values  of  Vshift  are  determined  by  the  size  and  shape  of 
the  region  I+II  and  can  be  easily  calculated.  The  slight  difference  in  region  I+II  centered  at  different 
rings  results  in  different  Vshift:  +4.398  eV  for  the  3-member  ring  and  +4.424  eV  for  the  6-member  ring. 
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The  value  of  Vshm  for  a-quartz  (-3.369  eV)  differs  dramatically  from  these  for  amorphous  silica  due  to 
a  very  different  shape  of  region  I+II  and  larger  density  of  quartz.  To  compare  values  of  Einc,  IP(H") 
and  EA(H+)  calculated  for  different  systems,  these  energies  have  to  be  corrected  by  VsiliftxAQ,  where 
AQ  is  the  change  of  the  total  charge  of  the  QM  cluster  in  a  corresponding  process. 

The  calculated  values  presented  in  Tables  1  and  2  are  corrected  for  Vshift.  The  Mott-Littleton  correction 
due  to  the  polarization  of  the  infinite  solid  outside  region  I  is  about  0.3  eV  and  is  not  included  into 
Tables  1  and  2.  The  frequencies  of  O-H  and  Si-H  stretching  modes  were  calculated  in  the  harmonic 
approximation  by  diagonalizing  the  dynamic  matrix.  The  latter  was  calculated  for  the  displacements  of 
the  Hydrogen  ion  only  with  the  electronic  structure  of  the  QM  cluster  being  relaxed  for  each  position 
of  the  H  ion. 


3.  Results  of  test  calculations 

As  was  noted  above,  one  of  the  aims  of  this  project  has  been  the  development  of  computational  tools 
for  analyzing  amorphous  structures.  The  construction  of  the  most  statistically  meaningful  amorphous 
structures  is  a  formidable  task  on  itself,  which  has  not  been  fully  resolved  so  far  [6].  The  structures 
discussed  in  this  chapter  are  obtained  using  the  best  recipes  suggested  in  the  literature  on  MD 
generating  of  amorphous  silica  structures  [16,17].  They  are  used  as  representative  systems  to 
demonstrate  the  applicability  of  the  developed  methods  to  study  defects  in  amorphous  silica. 

Several  glass  structures  generated  using  the  Molecular  Dynamics  technique  were  analyzed  in  terms 
of  distributions  of  inter-atomic  distances,  angles  between  neighboring  atoms,  and  on-site  electrostatic 
potential.  The  obtained  distributions  are  qualitatively  very  similar,  which  demonstrates  that  their 
statistical  properties  do  not  depend  significantly  on  the  parameters  of  MD  calculations  and  the  number 
of  atoms  in  the  simulation  box.  This  is  illustrated  in  Figs.  4  and  5  presenting  the  results  for  two 
structures  containing  192  and  648  atoms  in  the  simulation  box,  respectively.  The  shape  and  width  of 
both  distributions  are  very  similar  and  are  in  general  agreement  with  the  results  of  previous  theoretical 
modeling  [16,17]  and  with  the  experimental  data  [6]. 


Fig.  4.  Distribution  of  interatomic  distances  in  structures:  192  atoms  -  red  and  648  atoms  -  blue. 
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Using  the  PAS  code  we  have  also  analyzed  the  ring  statistics  in  five  MD  generated  structures.  Unlike 
a-quartz,  which  contains  only  6-  and  8-member  rings,  amorphous  silica  can  be  characterized  by  the 
ring  distribution,  which  includes  3-  to  8-member  rings.  Although  the  concentration  of  3-member  rings 
is  relatively  small,  they  are  different  from  larger  rings  and  therefore  may  have  some  special  properties. 
For  example,  as  was  already  mentioned  above,  they  are  thought  to  determine  the  tail  of  the  optical 
absorption  in  amorphous  silica  [13].  On  the  other  hand,  the  most  abundant  are  6-member  rings,  which 
are  also  present  in  all  crystalline  modifications  of  silica. 


O-Si-O  angles,  degrees 


Fig  5.  Distribution  of  Si-O-Si  and  O-Si-O  angles  for  structures:  192  atoms  -  red  and  648  atoms  -  blue. 

Future  defect  studies  will  require  more  detailed  analysis  of  amorphous  structures  for  specific 
properties.  For  example,  excitons  and  holes  may  localize  only  at  specific  sites.  To  find  these  specific 
sites  one  may  need  to  look  into  the  local  properties  such  as  correlations  of  the  ring  size  with 
interatomic  distances  and  angles  associated  with  atoms  in  the  ring.  An  example  of  such  a  correlation  is 
presented  in  Figs.  6  and  7.  Fig.  6  shows  the  distribution  of  distances  between  the  atoms  participating  in 
3-  and  6-member  rings.  These  distributions  were  calculated  in  the  following  way.  First,  information  on 
the  number  and  structure  of  rings  was  accumulated  as  discussed  in  the  previous  section.  Then  all 
atoms  associated  with  3 -member  rings  were  found  and  the  distribution  of  inter-atomic  distances  was 
calculated.  The  latter  step  of  the  calculation  was  repeated  for  6-member  rings.  Fig.  7  shows  the 
distributions  of  angles  calculated  in  a  similar  way. 
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Fig.  6.  Distribution  of  Si-0  and  Si-Si  inter-atomic  distances  for  3-member  rings  (red)  and  6-member  rings 
(blue)  within  a  648-atom  structure. 


Closer  examination  of  these  distributions  for  e.g.  3-member  rings  allows  us  to  deduce  the  following 
features: 


i)  Si-Si  distances  within  3-member  rings  are  at  the  lower  end  of  the  distribution,  i.e.  they  are  on 
average  about  0.10  -  0.15  A  shorter  than  the  Si-Si  distances  within  larger  rings.  It  can  also  be 
seen  that  the  maximum  of  the  Si-Si  distributions  shifts  to  larger  distances  for  larger,  6- 
member  rings.  However,  we  do  not  have  enough  statistics  for  more  quantitative  analysis. 

ii)  Si-O-Si  and  O-Si-O  angles  associated  with  3-member  rings  are  distinctly  smaller  than  the 
angles  associated  with  larger  rings.  It  is  possible  than  the  same  statement  also  applies,  but  to  a 
smaller  extent,  to  4-member  rings. 

iii)  There  is  no  obvious  correlation  between  the  ring  size  and  the  on-site  electrostatic  potential. 


O-Si-O  angles,  degrees 

Fig.  7.  Distribution  of  Si-O-Si  and  O-Si-O  angles  for  3-member  rings  (red)  and  6-member  rings  (blue)  within  a 
648-atom  structure. 
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To  test  our  embedded  cluster  method  for  defect  calculations  in  amorphous  silica  we  considered 
Hydrogen  defects  in  3-  and  6-member  rings  in  the  largest  648-atom  amorphous  structure.  The 
calculated  properties  of  H+  and  H~  centers  associated  with  a  3 -member  ring  and  a  6-member  ring  are 
summarized  in  Tables  1  and  2  respectively.  The  properties  of  these  centers  in  a-quartz  calculated 
using  the  same  method  are  also  included  for  comparison.  To  consider  the  effect  of  different  types  of 
local  environment  on  the  properties  of  H  centers,  we  report  a  range  of  their  characteristics  including 
formation  energies,  parameters  of  the  local  atomic  structure,  vertical  electron  affinities  and  ionization 
potentials,  harmonic  frequencies  of  H  ion  vibrations.  To  analyze  the  local  electronic  structure  of  the  H 
centers  we  also  report  atomic  charges,  which  demonstrate  strong  electron  density  redistribution.  The 
atomic  charges  are  given  for  Hydrogen  atom,  the  nearest  neighbor  O  or  Si  atoms  (in  the  case  of  H+  or 
H'  centers,  respectively),  and  for  one  of  the  second  nearest  neighbor  atoms.  The  values  of  vertical 
defect  ionization  energies  and  electron  affinities  were  calculated  in  two  ways:  with  frozen  shells 
(unrelaxed)  on  the  atoms  surrounding  the  QM  cluster  and  including  the  shell  relaxation  (relaxed).  The 
latter  simulates  the  effect  of  the  electronic  polarization  of  these  atoms  and  is  included  for  comparison. 


Fig.  8.  Local  atomic  strictures  of  the  H  and  H+  centres  in  both  a-quartz  and  a-Si02.  See  Tables  1  and  2  for 
comparison  of  the  geometrical  parameters. 

We  have  calculated  properties  of  the  H+  center  in  a-quartz  using  QM  clusters  Si207Si  6  and 
Si50i6Si*i2-  The  results  of  these  calculations  (see  Table  1)  suggest  that  the  dependence  of  the  results  on 
the  cluster  size  is,  indeed,  rather  small.  Thus  we  expect  that  the  difference  in  the  shape  and  size  of  QM 
clusters  used  in  calculating  the  defect  properties  in  a-Si02  will  not  bring  significant  uncertainty  to  our 
results  for  the  H+  centers  in  3-  and  6-  member  rings.  The  estimated  errors  are  about  0.2  eV  for 
formation  energies  and  about  0.1  eV  for  vertical  electron  affinities.  We  expect  that  similar  estimates 
will  be  true  for  the  H“  center  as  well,  but  more  studies  are  needed  to  confirm  this  assertion. 

Our  results  for  the  H+  center  suggest  that  the  defect  properties  depend  significantly  on  particular  site  in 
amorphous  structure.  Without  further  studies  we  could  not  establish  any  systematic  pattern  of  the 
dependence  of  defect  properties  on  the  ring  size.  Nor  can  we  state  which  of  the  two  rings  studied  can 
better  match  the  local  properties  of  a-quartz.  Indeed,  the  formation  energies  of  the  H+  center  in  a- 
quartz  and  the  3-member  ring  are  almost  identical,  while  the  6-member  ring  gives  the  much  lower 
formation  energy.  On  the  other  hand,  the  electron  affinities  of  the  H+  center  in  a-quartz  and  the  6- 
member  ring  are  very  close,  while  that  for  the  3-member  ring  is  almost  1  eV  higher.  Finally,  the  O-H 
stretch  vibrational  frequencies  are  similar  for  the  H+  center  in  the  3-  and  6-member  rings  and  by  about 
150  cm'1  higher  than  in  quartz. 

Comparison  of  the  O-H  and  Si-H  stretching  frequencies  for  different  H+  and  H~  centers  in  Tables  1 
and  2  shows  that  the  variations  of  Si-H  frequencies  are  much  smaller  than  the  variations  of  O-H 
frequencies.  Similarly,  the  differences  in  geometrical  parameters  (see  Fig.  8)  for  different  H“  centers 
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are  significantly  smaller  those  for  the  H+  centers.  This  can  be  rationalized  if  we  note  that  Si04 
tetrahedral  units  are  very  rigid.  Therefore,  a  defect  associated  with  such  a  unit,  H~  for  example,  is 
“protected”  from  the  rest  of  the  solid.  On  the  contrary,  the  H+  center  is  associated  with  an  Oxygen 
atom  connecting  two  Si04  units,  which  can  be  relatively  easily  rotated  with  respect  to  each  other. 
Therefore,  the  properties  on  H+  center  are  determined  by  relative  positions  of  at  least  two  tetrahedra. 


Table  1.  Properties  of  H+  centers  in  a-quartz  and  a-Si02. 


Property 

a-quartz 

3 -mem  ring 

6-mem  ring 

Si207Si  6 

Si50i6Si  12 

SLOgSi  6 

Si60i8Si  12 

^incorporation) 

-9.374 

-9.588 

-9.521 

-14.467 

Local  geometry, 

d(O-H) 

0.985 

0.983 

0.972 

0.967 

A  /  deg 

p(Si-O-Si) 

134.6 

135.2 

130.2 

122.1 

Vertical  EA, 

unrelaxed 

2.474 

2.367 

3.037 

2.249 

eV 

relaxed 

3.048 

2.891 

3.606 

2.727 

co  (H+),  cm'1 

O-H  stretching 

3530.7 

3602.1 

3740.9 

3788.6 

Atomic  charges 

H 

0.645 

0.636 

0.617 

0.614 

(ground  state). 

0 

-1.189 

-1.195 

-1.177 

-1.192 

M 

Si 

2.845 

2.834 

2.803 

2.835 

Atomic  charges 

H 

0.609 

0.601 

0.603 

0.594 

(after  e-  trapping), 

0 

-1.174 

-1.184 

-1.166 

-1.181 

M 

Si 

2.364 

2.359 

2.285 

2.332 

Table  2.  Properties  of  H  centers  in  a-quartz  and  a-Si02. 


Property 

a-quartz 

3 -mem  ring 

6-mem  ring 

Si50i6Si  12 

S^OgSi  6 

SieOjgSi  12 

Eincorporation) 

3.324 

4.931 

9.247 

Local  geometry, 

d(O-H) 

1.455 

1.449 

1.442 

A  /  deg 

P(Oax-Si-Oax) 

169.8 

172.8 

168.7 

P(Oeq-Si-Oep) 

106.2 

107.6 

107.0 

Vertical  IP, 

unrelaxed 

7.189 

6.758 

7.129 

eV 

relaxed 

6.741 

6.299 

6.772 

co  (H“),  cm'1 

Si-H  stretching 

2336.6 

2347.8 

2389.3 

Atomic  charges 

H 

-0.255 

-0.253 

-0.286 

(ground  state), 

Si 

2.340 

2.344 

2.387 

M . . 

O 

-1.315 

-1.320 

-1.308 

Atomic  charges 

H 

-0.007 

-0.151 

-0.070 

(ionised  state), 

Si 

2.256 

2.277 

2.332 

_ tel _ 

O 

-0.772 

-0.586 

-0.708 

We  would  like  to  point  out  that  the  H  centers  discussed  above  do  not  only  differ  by  their  environment 
Different  QM  clusters,  both  in  terms  of  their  topology  and  size,  were  used  in  our  calculation.  The  latter 
must  contribute  to  some  extent  to  the  spread  of  values  presented  in  Tables  1  and  2.  However,  from  our 
experience  with  the  neutral  Oxygen  vacancy  in  a-quartz,  the  dependence  of  the  formation  energies  on 
the  cluster  size  is  quite  small  [11].  Therefore,  preliminary  we  can  conclude  that  variations  in  properties 
of  localized  defects,  such  as  H+  and  H'  centers  at  different  sites  of  Si02  polymorphs  and  in  a-Si02,  are 
due  to  differences  in  their  environment  at  ranges  II  and  III.  The  local  atomic  structure  and  charge 
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density  of  these  defects  are  essentially  independent  on  the  details  of  their  environment.  Finally,  we 
note  that  an  ^-member  ring  can  potentially  accommodate  n  different  H+  centers  and  up  to  6n  different 
H"  centers.  Therefore,  further  calculations  are  needed  to  understand  better  whether  properties  of 
localized  defects  such  as  FT  and  H“  centers  correlate  with  the  ring  size  and,  more  generally,  with  their 
short-,  medium-,  and/or  long-range  environment  (ranges  I,  II,  and  III  respectively). 

4.  Summary  and  conclusions 

In  this  project  we  developed  computational  tools  for  analysis  of  amorphous  structures  and  extended 
the  embedded  cluster  approach  to  study  localized  defects  in  amorphous  Si02.  The  PAS  code  allows 
one  to  fully  characterize  amorphous  structures  in  terms  of  distributions  of  bond-lengths,  angles,  rings 
and  other  structural  parameters.  It  also  allows  one  to  correlate  local  geometric  parameters  of 
amorphous  structures  with  electrostatic  potential  and  search  for  particular  structural  features.  The 
embedded  cluster  method,  in  its  turn,  allows  one  to  perform  accurate  quantum  mechanical  calculations 
of  selected  defect  sites  and  their  local  environment  and  to  account  classically  for  the  short-  and  long- 
range  response  of  the  rest  of  the  amorphous  network.  This  method  is  self-consistent  in  a  sense  that  a 
QM  cluster  is  built  as  a  fragment  of  an  amorphous  structure  and  remains  coupled  to  the  amorphous 
network  throughout  the  calculation.  Moreover,  the  classical  environment  of  the  QM  cluster  is 
described  using  the  same  set  of  inter-atomic  potentials,  which  was  used  to  generate  the  amorphous 
structure  in  Molecular  Dynamics  simulations. 

We  have  built  several  amorphous  structures  using  classical  Molecular  Dynamics  and  used  the  PAS 
code  to  analyze  their  properties  and  to  search  for  correlations  between  these  properties.  The  results  of 
this  analysis  demonstrate  that  our  structures  are  statistically  equivalent  to  those,  obtained  in  previous 
calculations.  For  the  first  time  we  applied  the  embedded  cluster  method  to  study  defects  centers  in 
amorphous  Si02.  The  H+  and  H~  centers  associated  with  small  (3-member)  and  large  (6-member) 
rings  were  considered  and  their  properties  were  calculated.  The  results  of  suggest  that  there  are 
considerable  site  variations  of  the  FT  and  H“  defect  properties.  These  variations  are  larger  than  those 
induced  by  the  size  and  shape  of  QM  clusters  used  in  the  calculation  and  therefore  can  be  studied 
using  the  embedded  cluster  method  developed  in  this  project. 
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